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A BRIEF DESCRIPTION OF THE JAMESON- CAUGHEY 
NYU TRANSONIC SWEPT-WING COMPUTER 
PROGRAM - FLO 22 

Antony Jameson (NYU) 

David A. Caughey (Cornell) 

Perry A. Newman (NASA) 

Ruby M. Davis (NASA) 

Langley Research Center 
PREFACE 

This document was prepared by the third and fourth authors 
to serve as an interim user guide for the Jameson- Caughey NYU 
Transonic Swept-Wing Computer Program - FLO 22. This information 
is pertinent to the version provided to NASA LaRC in the spring 
of 1976. The first and second authors are in the process of 
preparing more extensive documentation on a later version of 
program FLO 22. This document is intended to meet several 
immediate demands for use-e-f this computer code in NASA-sponsored 
programs. 

The first two sections entitled "Calculation of the Flow 
Past a Swept Wing" and "Some Results of Swept Wing Calculations" 
are sections 5.7 and 5.8 of reference 1 wich very minor editing. 
The last section entitled "Input Description" is an update of 
reference 2 (an input description for Jameson's yawed wing 
program - FLO 17) which includes the swept-wing input parameters. 


SUMMARY 


Prof. Antony Jameson (NYU) and Prof. Dave Caughey (Cornell) 
have developed this computer program for analyzing inviscid, iser- 
tropie, transonic flow past 3-D swept-wing configurations. This 
work was done at the Courant Institute of Mathematical Sciences, 
New York University, under NASA Grants NGR-33-016-167 and 
NGR-33-016-201 . Some basic aspects of the program are: The 

free-stream Mach number is restricted only by the isentropic 
assumption. Weak shock waves are automatically located where 
ever they occur in the flow. The finite-difference form of the 
full equation for the velocity potential is solved by the method 
of relaxation, after the flow exterior to the airfoil is mapped 
to the upper half plane. The mapping procedure allows exact satis 
faction of the boundary conditions and use of supersonic free 
stream velocities. The finite difference operator is "locally 
rotated" in supersonic flow regions so as to properly account 
for the domain of dependence. The relaxation algorithm has been 
stabilized using criteria from a time-like analogy. The brief 
description contained in this document should enable one to 
use the program until a formal user's manual is available. 


CALCULATION OF THE FLOW PAST A SWEPT WING 
It is desired to solve the three-dimensional potential flow 
equation which can be written in quasilinear form as 

(a 2 -u 2 )<O xx + (a 2 -v 2 )<f> yy + (a 2 -w 2 )<i> zz 

- 2uvtf xy “ 2vw4> yz " 2uw *xz “ 0 ' (1) 

where u, v and w are the velocity components and a is the 
local speed of sound. The singularity at infinity in the velocity 
potential is removed by introducing a reduced potential 

G ■ <f> - x cos a - y sin a (2) 

where a is the angle of attack. 
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In the case of a lifting flow the velocity potential 
is discontinuous across the vortex sheet trailing behind 
the wing. Roll up of the vortex sheet will be ignored; 
the conditions to be satisfied at the surface in which 
the vortex sheet is assumed to lie are that the jump 

T in the potential is constant along lines parallel to 

the free stream, and that the normal component of velocity 

is continuous through the sheet. At infinity the i.-ow is 
undisturbed except in the Trefftz plane far downstream- 
where there will be a two dimensional flow induced by the 
vortex sheet. 

The construction of a satisfactory curvilinear 
coordinate system to suit the geometry of the configuration 
is one of the most difficult aspects of the three dimensional 
problem. Here nonorthogonal coordinates will be generated 
by a sequence of elementary transformations. First parabolic 
coordinates are introduced in planes containing the wing section 
by the square root transformation 


1/2 

x i + iY i B C x “ x 0 («> + My-y 0 (z))] , z x * z , 

where z is the spanwise coordinate, and Xq and define 
a singular line of the coordinate system located just inside 
the leading edge (see Figures 1 and 2). The effect of this 
transformation is to unwrap the wing to form a shallow bump 




Iff 


n 




Then a 


Y x * S(X 1 ,Z 1 ) 

shearing transformation is used 


(d) 


X - X x , Y = Y x - S(X 1 ,Z 1 ) , Z = Z 


(5) 


to map the wing surface to a coordinate surface. Finally, 
in order to obtain a finite computational domain X, Y and Z 
are replaced by stretched coordinates X, Y and Z. 

The stretching used in the present computer program is to 
set X « X in an inner domain -X < X < X , and to set 

lu ““ XU 

* - *m + < S ‘ V f - (l^) 2 } «) 

rn 

when X > X ffl , with a corresponding formula when X <-X jn , 
so that X = ±® -when X * ±1. Typically the parameter 
a has the value 1/2. Similar stretchings are used for 
Y and Z. 

The vortex sheet is assumed to coincide v?ith the cut 
behind the singular line which is opened up by the square 
root transformation (Eq. (3)). Thus a jump T is introduced 
in the potential between corresponding points representing 
the two sides of the vortex sheet. A complication is 
caused by the continuation of the cut beyond the wing. 

Points on the two sides of the cut must be identified as 
the same point in the physical space. Also a special form 
of the equations must be used at points lying on the singular 
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line beyond the wing. At these points the equation to be 
satisfied reduces to the two dimensional Laplace equation 


in- the and Y^ coordinates. An advantage of the square 
root transformation (Lq. (,!•»)) is that i i collapses the height 
of the disturbance due to the vortex sheet to zero in the 
parabolic coordinate system at points far downstream, 
where approaches infinity, with the result that 

the far field boundary condition is simply 

G = 0 (7) 

The final equation for the reduced potential G contains 
numerous terms. In order to construct a rotated difference 
scheme with a proper upwind bias at supersonic points it is 
only necessary, however, to consider the principal part, 
consisting of the terms containing the second derivatives of G. 

Suppose that equation (1) is written in the canonical form 

(a 2 - q 2 M> ss + a 2 (A<j> ~ 4> ss ) - 0 (8) 

3 2 a 2 a 2 

where A is the Laplacian operator — + — *• + — and 

ax'* 3y z 3z^ 

d> is the streamwise second derivative in a Cartesian 
T ss 

coordinate system locally aligned with the flow 

♦s, - -h (» 2 *xx + v %y + w \z + Juv *xy + 2vw<, yz + 2u "’^z) • (9) 

Then at supersonic points upwind difference formulas are used 
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for all second derivatives of G arising fror. the transforr.a- 


tion of $ 


ss 


into the curvilinear coordinate system, and 


central difference formulas for all second derivatives of G 
arising from the- transformation of - <? gs * The formulation 
in terms of the Laplacian avoids the need to determine 
explicitly a pair of local coordinate directions normal to 
the stream direction. 


The difference equations are solved by relaxation, 
with care taken to make sure that at supersonic points the 
equivalent time dependent equation is compatible with the 
steady state equation. If m and n are coordinates in 
a plane normal to the streamwise direction s, and M is 
the local Mach number q/a, the equivalent time dependent 
equation can be written in the form 


<K 2 -1> V W 2o lW 2o 2*nt + 2a 3^nt + ’ 1 

where the coefficients , a 2 and q 3 depend on the split 
between new and old values used in the relaxation scheme. 

To make sure that this is a wave equation with s as the time- 
like direction, an analysis indicates that the difference 
formulas should be organized so that 


“1 > 


/(K Z -1) ( q 2 +ci 3^ 


(11 


8 


Also the balanced coefficient rule should still be 
applied in the supersonic zone, corresponding to a 
2 ero coefficient of <j> t in equation (10). It is 
convenient to solve the equations for the correction to the 
potential simultaneously along lines, corresponding to 
a point relaxation process in two dimensions. Any of the 
coordinate lines can be used for this purpose, the choice 
being guided by the need to avoid advancing through the 
supersonic zone in a direction opposed to the flow. In 
practice it has been found convenient to divide each X, Y 
plane into three strips, and to march towards the wing 
surface in each central strip, updating horizontal lines, 
and then outwards in the left-hand and right-hand strips, 
updating vertical lines. (See Figure 3.) 


MARCHING DIRECTION 


HORIZONTAL LINE RELAXATION 

FIGURE 3. MARCHING DIRECTIONS OF RELAXATION SCHEME 
FOR SWEPT KING CALCULATION 
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SCLIE RESULTS OF SWEPT WING CALCULATIONS 




Some results of swept wing calculations are presented 
in Figures 4 ancl 3. These were calculated on a grid 
with 144 cells in the chordvise X direction, 24 cells in 
the normal Y direction and 32 cells in the spanwise Z 
direction, calling for the solution of the difference equations 
at 109824 mesh points. In each case the interpolated 
result of a preliminary calculation on a 72* 12x16 grid 
was used to provide the starting guess: 200 cycles were 
used on the coarse grid followed by 100 cycles on the fine 
grid. Such a calculation requires about 80 minutes on a 
CDC 6600 (which should be reduced to about 20 minutes on a 
CDC 7600) . 

Figure 4 shows the result of a calculation for a 

rather simple wing tested by ON ERA, for which experimental 
data has been published. It can be seen from Figures 4(c) 
and 4(d) that the agreement with the experimental 
data is quite good, despite the fact that no attempt was 
made to allow for viscous effects. As in the case of the 
two dimensional calculations, the nonconservative difference 
scheme introduces a source distribution over the shock 
surfaces, causing a displacement of the streamlines and a 
forward shift in the location of the shock waves. Apparently 
this partially compensates for the absence of a correction 
for the displacement effect of the boundary layer. 


I I 
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To illustrate the geometric complexity afL the 
configurations which can be treated by the program, Figure 
5 shows the results of a calculation for a wing designed 
and tested by the Douglas Aircraft Company. The wing 
is a typical design for a long range transport aircraft, with a 
sweepback of 35 degrees at the leading edge, and a substantial 
change in the-section between the root and tip. The. test 
was of a wing-body combination. In the calculation the v/ing 
was extended to the plane of symmetry at the., fuselage center - 
line. The calculation shows two shock waves over the 
inboard part of the wing. The forward shock wave originates 
from the leading edge at the wing root, and the aft shock 
wave is roughly normal to the free stream. The two shock 
waves merge at about the 1/4 span point, forming a triangular 
shock pattern over the upper surface of the wing. The 
coalescence of the shock waves can be traced in Figures 
5(c) - 5(g), in which the pressure distributions at 

a sequence of span stations over., the- inboard part of the 
wing are plotted separately. The convergence history of 
this calculation, measured by the largest residual, is 
shown in Figure 6. 




12 



VIEW OF WING 


ONERA WIN 3 He 

L.E 

MACH .923 

YAW 

L/D -.00 

CL 


0.000 

-.0000 


ALPHA O.C 
CD .02M6 


FIGURE 4(A) 


ry 


PER SURFACE PRESSURE 


LONER SURFACE PRESSURE 


OPEr.P, K'lliG M6 

L.E. 

SICFi' 30 

DE3 ASPECT RATIO 

ARCH 

«9r?3 

YflW 

o.cvo 

ALPHA 0.000 

L/D 

-.00 

CL 

- .0030 

CD .0246 


FkillHI 'HID 





» 



OPERA WHIG M6 L.E, SWEEP 30 DEG ASPECT RATIO 

MACH .923 YAW 0.000 ALPHA 0 

2 ,80 CL - .0000 CD 

• Tin-'rv LXPIIUHEHT 

MGIIHI 4(D) 








VIEW OF WING 
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INPUT DESCRIPTION 


Read Number 

Order, Cards . 

1 1 


2 1 

3 1 


Lefccription and Comments 


title. - Descriptive title of case or sequence; 

Format (8A10) . . 

Appears on Varia in Plots and beginning 

of output. 


DESC. - Description for card in Read Order 3 
Format (8A10). 


FNX, FNY, 


FNZ, FPLOT , XSCALE, PSCALE, FCONT 
Format ( 8E10 . T ) 

Note: A number cf quantifies are read in 

as floating-point numbers and converted 
tc integers -within the program. 


FNX. 


m. 


m. 


- Number of computational grid points in 
"cbordvise direction" from downstream 
infinity, around the leading edge and 
bach to downstream infinity on cof rses t^ 
mesh. Maximum is 96 (192 with nc grid 
halving). 


Number of computational grid points in 
"nor m* 1 direction" from airfoil surface 
to infinity on coarsest mesh . Maximum is 
12 ( 24 with no grid halving). 

Number of computational grid points in 
"soanwise direction" from infinity, 
across the wing span and to infinity on 
coarsest mesh. Maximum is 1 6 (52 with no 


FPLOT. - Plot trigger. Selects type of plot for 
chordwise surface pressure coefficients. 
FPLOT » 0. Printer plots, one at each 
soanwise grid plane section with C? 
versus the computational grid chordwise 
variable. 


25 


Read Number Descript icr. as A Comments, 

Ord er FPLOT * 1* Varian plots (from THREED). 

These are superimposed plots, vith all 
span sections shown or. two figures, an 
upper surface and a lower surface plot 
of CP versus physical space chord- 
wise variable. 

FPLOT ■ 2. Varian plots (from THREED) 
as above plus section plots (from GRAPH.)* 
These latter plots-, one per section, 
give upper sued lower surface C? versus 
physical space chordvise variable. 


Defaults to zero 


XSCALE. - Scale of abcissa in pressure plots at each 
wing station: XSCALE > 0. » Length of 

section airfoil chord (in inches); 

XSCALE * 0. , 5-inch section chord; 

XSCALE < 0. * Length of maximum section 
airfoil chord will be (XSCALE! inches on 
plot. All other section plots will be 
scaled according to wing planform. 


k 

5 



P SC ALE. - 


FCONT. - 


1 

DESC. - 

1 card 

FIT, COVO, 

for each 


computational 


grid. Maximum 

FIT. 

essentially 2, 


dimensioned J. 

COVO. - 


Scale of ordinate in pressure plots, C^ 
on plot: PSCALE ■ 0.# -.4 per inch. 

PSCALE + 0., - |PSCALE| per inch. 

Program starting/continuing trigger. 

FCONT * 0. t Calculation begins at itera- 
tion zero. FCONT * 1. , Calculation is to 
be continued from a previous one. 


Description for card in Read Order 5 
Format (8A10). 

P10, P20, PJ0, 3ETAC, STRIP0, FHALF 
Format (8E10.7) 

Maximum number of iterations on 
this gird, called MIT in program. 

Convergence criterion on the maximum 
change in reduced velocity potential 
(G) from one iteration cycle to the 
next on this prid. 
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Read Number 

Order Cards 


Description and Commentc 


P10. - Subsonic point relaxation factor on 

this grid; must be <2.. Typically 
1.6 on coarse grid. 

F20. - Supersonic point relaxation factor; 

must be & 1.. Should use 1.0 for 
stability on all grids. 


6 1 

7 1 




P30. 

- Circulation relaxation factor, May 
be > 1.0 . 

beta: . 

- Stabilization factor used at super** 
sonic points in finite difference, 
operator if 3ETA0 > 0. Most needed 
when *1., many cases operate satis- 

factorily with 3ETAC = 0. Convergence 
is slowed but stability enhanced when 
3ETAO0. 

STRIPC. - 

Line relaxation control. Computational 
X-Y planes are relaxed by horizontal 
lines (YSWEEP) in central strip, verti- 
cal lines (XSWEEP) in outer strips. 
STRIPO specifies the fraction of compu- 
tations 1 plane included in central 
strip: 0. £ STRIPO < l., where STRIPO-1. 

all horizontal line relaxation. 

FHALF. - 

Grid halving trigger. 

|FHALF|> 1. read another card (Read 
Order 5 fonnat) containing computation- 
al parameters to be used on grid with 
mesh size halved in all directions. 
|FHALF|< 1_ must appear on finest grid 
card (last one read).- Calculation pro- 
ceeds automatically through the sequence 
of computational grids. 

DESC. 

Description for card in Read Order 7 
Format (8A10). 

FMACH, YA, 

AL, CDO 

Format (8EIO.7) 

FMACH. - 

Free stream Mach number. 

YA. 

Yaw angle (in degrees). 

AL. 

Angle of attack (in degrees) measured 
in plane containing freeatream direction. 

CDO. 

Drag coefficient due to skin friction 
(CD FRICTION on output). This input 
number is added tc the drag coefficient 
obtained by integrating the surface 
pressures (CD FCRK on output). 
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Read Humber Description and Comments 

Order Cards 

Read Orders 8 through 19 are used to specify the wing geometry (in physical 
space, of course). One can define the wing at up to 11 span stations. A 
set of airfoil coordinates must be read in at the first station. It need not 
be read in at other stations, if one is changing -only combinations of the 
following three airfoil section parameters: chord., thickness ratio or 

angle of attack (twist) The wing -shape at intermediate span positions 
(i.e., the computational grid planes for example) is obtained by linear. - 
interpolation in the spanwise direction in the physical space. 


Read Orders 8 and 9 are read only once: 10 and 11 are read- FNC (see 9) times.; 

12 through IT (19 if non- symmetric airfoil section) must be read at first 
section and may be required at other sections, .depending on the wing geometry. 

8 1 DESC. - Description for card in Read Order 9 

Format (8A10), 

9 1 ZSYM, FNC, SWEEP 1, SWEEP 2, SWEEP, DIKED 1, 

DIKED 2, DIHED 


ZSYM. - Wing planform symmetry trigger: 

ZSYM ** 0 * yawed wing 
ZSYM = 1 , swept wing 

FNC.- - Must be ~ 3 . The leading edge of the wing 
in physical space is fit with a cubic 
spline. Data at three span stations are 
required (as minimum) as wall as the six 
angles which follow. If the wing leading 
edge has a slope discontinuity, three 
stations should be used fairly close to 
it. 

SWEEP 1.- Sweep angle of wing leading edge at root 
section (in degrees) . 

SWEEP 2.- Sweep angle of wing leading edge at tip 
• section (in degrees). 

SWEEP. - Sweep angle of spanwise grid lines at 00 
(off tip of wing) (in degrees). 


DIHED 1^«-Dihedral angle of wing leading edge at root 
section (in degrees). 

DIHED 2. - Dihedral angle of wing leading edge at tip 
section (in degrees). 


- Dihedral angle of spanwise grid lines at 
(off tip of wing) (in degrees). 
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DIHED. 


00 


Read 

Order 


Number 

Cards 


Description and Comwont. 


10 

1 

DSSC. 

11 

1 

ZS (K) 



’ ZS(K) 


XL. 

YL. 

CHORD. - 
THICK. - 

AL. 

FSEC. - 


Description for cards in ?.= ad. -Crier il 
Format 48A10). 

YL, CHORD, THICK, AL, FSEC 
Format- (8EI0. 7) 

Spanvise. coordinate of the wing section 
being specified— It. is in the same 
uxdts~as CHORD. These stations are 
ordered from tip-to-tip, ia ascending 
algebraic order, of ZS(K) for yawed wing 
and root-to-tip for swept wing. 

X coordinate of section leading edge in 
physical space (controls sweep). 

Y coordinate of section leading edge in 
physical space (controls dihedral). 

Section chord length. The chord of 
the airfoil coordinates to be read in 
(or already read in at the prior 
station) will be scaled to this 
value. 

Section thickness ratio relative to 
that of the airfoil coordinates to be 
read in (or already read in- at the 

prior station). Note, this is a ratio 
of thickness/chord ratios. The thick- 
ness of the airfoil coordinates will 
be scaled with this Value- 

Section angle of attack or twist 

(in degrees). Airfoil coordinates will 

be rotated through this angle about LE. 

Section airfoil coordinate trigger. 

FSEC * 0. Do not read airfoil coor- 
dinates. Last set of airfoil coordin- 
ates read will be used at this section. 
They may be scaled by any combination 
of' CHORD, THICK, or AL read above. 
Skip Read Orders 12 through lj? for 
this section. 


Bead 

Order 

Number 

Cards 


Description and Comments 




rSEC = 1. Read a new set cf airfcil 
cccrdinates which will be used at this 
station ana perhaps at ether stations. - 
They may be scaled by any combination 
cf CHORD, THICK, or A I read above 

for this section^. 

At first station (K = 1) FSEC is 
ignored; one must supply Read Orders 
12 through 17. 

12 

1 

D2SC. 

Description for cards in Read Order 13 
Fermat (8A10). 

13 

1 

FSYM, FNU, 

FNL 

Format (8E10.7) 



FSYM. 

Airfoil symmetry trigger. 




FSYM >_ 1. Symmetric airfoil. 

Read in only upper surface airfoil 
coordinates, ordered leading edge 
to— trailing edge. 




FSYM <1. Non-symmetric airfoil. 
Read in upper and lover surface 
airfoil coordinates, respectively, 
each set ordered leading edge to 
trailing edge. Note thet leading- 

edge points are 
included in-both sets— 



FNU. 

Number of coordinates read in for 
upper surface of airfoil. 



FKL. 

Number of coordinates read in for 
lower surface of airfoil. 

Ik 

1 . 

DESC. - 

Description for cards in Read Order 13 
Format (8A10). 


15 1 TRL, SLT, XSriS, YSING 

Format (8E10.7) 

TRL. - Included angle of trailing edge of air- 
foil (in degrees). For blunt TE it is 
upper slope angle minus lower slope angl£. 

SLT. - Slope of airfoil mean camber line at 
trailing edge. 


P.ead 

Order 

Number 

Caras 


Description and Comments 



>;s i ;:g . 

X coordinate of the origin of the 
mapping referenced to the airfoil 
leading edge. Recommend approximately 
X(LE) + 1/2 leading edge radius where 
R^j. is in the same units as XP(I) read 

below. 




- Y coordinate of the origin of the 
mapping referenced to the airfoil, 
leading edge. Recommend approximately 
Y(LE). 

„ * 
-c 

1 

DESC. 

- Description for cards in Read Order 17 
Format (6A10). 

1 - 

FNU 

X? I}, 

YP(I) 

Fermat (oElC.?/ 

- X coordinate of airfoil upper surface, 
ordered leading edge tc trailing edge. 




- Y coordinate of airfoil upper surface, 
ordered leading edge tc trailing edge. 
Note that there is only one pair of 
coordinates per card. 


If airfoil section is not symmetric ( FSYM<L) the airfoil lower surface 
coordinates must he read here. For symmetric airfoil (FSYK>1.), skip the 
tvc Head Orders 18 and 19 . 


1 DESC. - Description for cards in Read Order 19 

Format (8A10). 

FNL VAL, DUM - Format (8E10.7) 

VAL. - X coordinate of airfoil lover surface, 
ordered leading edge to trailing edge. 

DUM. - Y coordinate of airfoil lover surface, 
ordered leading edge to trailing edge. 
Note that there is only one pair of 
coordinates per card. 

Read Orders 10 through 19 complete the input for one span station. As 
indicated above Read Order 8, at least Read Orders 10 and 11 must be repeated 
for the remaining FNC-1 sections when FMC i 3. 


18 

19 
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Description and Comments 


Read Number 

Order Cards 


The program terminates by reading the first three Read Orders with FNX<1.} 
that is,. last three cards for a normal stop should be: 

1 1 TITLE . - End of Calculation 


2 1 


DESC. 


Description for card in Read Order 3 


3 1 


0 . 
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